function [out] = deskew(in)
[m,n] = size(in);
[x,y] = meshgrid(1:n,1:m);
w = in(:);
w = w/sum(w);
ux = w'*x(:);
uy = w'*y(:);
sx = x(:)-ux;
sy = y(:)-uy;
xx = w'*(sx.*sx);
yy = w'*(sy.*sy);
xy = w'*(sx.*sy);
[V,D] = eig([xx,xy;xy,yy]);
%dd = diag(D)
%[~,I] = sort(dd,'descend');
%R = V(I,:);
V = V.*repmat(sign(V(2,:)),2,1);
if (V(1,2)>V(2,2)) 
    V = V*[0,1;-sign(V(1,1)),0];
else
    V = V*[sign(V(1,1)),0;0,1];
end
R = V;
nx = x(:)-(n+1)/2;
ny = y(:)-(m+1)/2;
rv = R*[nx,ny]';
rx = rv(1,:)+ux;
ry = rv(2,:)+uy;
out = interp2(x,y,in,rx,ry,'near');
out = reshape(out,m,n);
end